/*
/*smog_vehicle_emissions.do: This program produces a figure showing the average 
emissions levels of various model years, and an appendix table relating emissions to inspection
type (particularly initial inspections of out-of-state vehicles*/

/* Top portion commented out because it requires data not available in the archive*/

clear all

capture log close
set more off

cd "/home/work/projects/smog_check/Smog Check Data"

use quarter county station_id svfr_ratio using asmair_star

duplicates drop
tempfile temp
save `temp'

//local debug if substr(vin,-3,2)=="11"

use fpr_precollapse `debug'

joinby station_id quarter using `temp'

gen asm = test_type == "A"

gen directed = inlist(insp_reasn,"D","S","P")

gen age = year -modelyr

drop if age <0

replace timesince = timesince/365

replace insp_reasn = substr(insp_reasn,-1,1)

gen biennial = (insp_reasn == "P" | insp_reasn == "B" | insp_reasn == "D" | insp_reasn == "S" | insp_reasn == "1")

gen coo = (insp_reasn =="C"| insp_reasn=="2")

gen init_reg= (insp_reasn== "I" | insp_reasn == "3")

drop if biennial + coo + init_reg == 0

gen othergen = modelyr > 1985


bysort vin (date teststart): gen cycle = sum(initial_inspection)

sort vin cycle date teststart

by vin cycle: gen re_pass = fail == 0 & sum(fail)>0  & _n == _N

xtset make
local controls  ib10.age ib152.quarter estmiles timesince asm  i.county

eststo hc: xtreg hc (1.coo 1.init_reg)#othergen `controls' if initial_inspection & year >= 1998 & modelyr>=1976, fe vce(robust)
eststo no: xtreg no (1.coo 1.init_reg)#othergen `controls' if initial_inspection & year >= 1998 & modelyr>=1976, fe vce(robust)
eststo co: xtreg co (1.coo 1.init_reg)#othergen `controls' if initial_inspection & year >= 1998 & modelyr>=1976, fe vce(robust)
eststo fail: xtreg fail (1.coo 1.init_reg)#othergen `controls' if initial_inspection & year >= 1998 & modelyr>=1976, fe vce(robust)

#delimit ;

esttab hc no co fail using fpr_initial.tex, replace se booktabs
mtitles("Ln(HC)" "Ln(\nox)" "Ln(CO)" "P(Fail)") 
indicate("Vehicle Age FE = *.age" "Calendar Quarter FE = *.quarter" "County FE = *.county" "VIN Prefix FE = _cons")
transform(estmiles @*100000 @*100000 )
coeflabels(1.coo#0.othergen "Change of Ownership \\ \quad\quad Model Years 1976-1985"
1.coo#1.othergen "\quad\quad Model Years 1985+" 1.init_reg#0.othergen "Initial Registration \\ \quad\quad Model Years 1976-1985"
1.init_reg#1.othergen "\quad\quad Model Years 1985+" estmiles "Odometer (00000s)"
timesince "Years Since Last Inspection (Years)" asm "ASM Test" directed "Directed Inspection" )
title(Emissions levels by Inspection Reason\label{tab:initial})
;
#delimit cr

foreach pollutant in co no hc{
	replace `pollutant' = exp(`pollutant')
}
replace co = co*10000


preserve
gcollapse co no hc if test_type=="A", by(modelyr fail) 
save modelyr_emissions.dta, replace*/
use modelyr_emissions.dta, replace

#delimit ;

twoway connect no modelyr if fail==1, yaxis(1) || connect no modelyr if fail==0, yaxis(1) lpattern(dash) msymbol(S) ||  if modelyr>=1975 ,
graphregion(color(white)) ytitle("Combined NOx ppm", axis(1)) xlabel(1975(5)2010)
name(no, replace) xtitle(Model Year)
legend(label(1 "Failed Inspections") label(2 "Passed Inspections"))
;

twoway connect co modelyr if fail == 1 || connect co modelyr if fail == 0,  msymbol(S) || if modelyr>=1975 & modelyr<2009 ||,
graphregion(color(white))  ytitle("Combined CO ppm") name(co, replace) xlabel(1975(5)2010)
legend(label(1 "Failed Inspections") label(2 "Passed Inspections"))  xtitle(Model Year);

graph combine co no, graphregion(color(white)) xsize(10);

graph export modelyr_emissions.pdf, replace;
/*
restore;



preserve;

gcollapse fail, by(modelyr year) fast;

save failrate.dta, replace*/
use failrate.dta, replace

#delimit ;
twoway connect fail modelyr if year == 2002 ||
connect fail modelyr if year == 2009 || if modelyr >=1975 & modelyr < 2009,
 xlabel(1975(5)2010) xtitle(Model Year)
 legend(label(1 "Inspections in 2002") label(2 "Inspections in 2009"))
graphregion(color(white))  ytitle("I/M Failure Rate") ylabel(0(.05).3)
;

graph export my_failrate.pdf, replace;

instanttex using mygraph.tex, replace fig(modelyr_emissions);
